

<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
  <meta charset="utf-8">
  
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  
  <title>Step 2: Classical MD simulation of fracture in Si &mdash; quippy 7ea04c5+ documentation</title>
  

  
  
    <link rel="shortcut icon" href="../_static/favicon.ico"/>
  
  
  

  
  <script type="text/javascript" src="../_static/js/modernizr.min.js"></script>
  
    
      <script type="text/javascript" id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
        <script type="text/javascript" src="../_static/jquery.js"></script>
        <script type="text/javascript" src="../_static/underscore.js"></script>
        <script type="text/javascript" src="../_static/doctools.js"></script>
        <script type="text/javascript" src="../_static/language_data.js"></script>
        <script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
        <script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"processClass": "math|output_area", "processEscapes": true, "ignoreClass": "document", "inlineMath": [["$", "$"], ["\\(", "\\)"]]}})</script>
    
    <script type="text/javascript" src="../_static/js/theme.js"></script>

    

  
  <link rel="stylesheet" href="../_static/css/theme.css" type="text/css" />
  <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
  <link rel="stylesheet" href="../_static/theme_overrides.css" type="text/css" />
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
    <link rel="next" title="Step 3: LOTF hybrid MD simulation of fracture in Si" href="adaptive-qmmm-step3.html" />
    <link rel="prev" title="Step 1: Setup of the Silicon model system" href="adaptive-qmmm-step1.html" /> 
</head>

<body class="wy-body-for-nav">

   
  <div class="wy-grid-for-nav">
    
    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
      <div class="wy-side-scroll">
        <div class="wy-side-nav-search" >
          

          
            <a href="../index.html" class="icon icon-home"> quippy
          

          
            
            <img src="../_static/hybrid.png" class="logo" alt="Logo"/>
          
          </a>

          
            
            
              <div class="version">
                7ea04c5+
              </div>
            
          

          
<div role="search">
  <form id="rtd-search-form" class="wy-form" action="../search.html" method="get">
    <input type="text" name="q" placeholder="Search docs" />
    <input type="hidden" name="check_keywords" value="yes" />
    <input type="hidden" name="area" value="default" />
  </form>
</div>

          
        </div>

        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
          
    
            
            
              
            
            
              <ul class="current">
<li class="toctree-l1"><a class="reference internal" href="../intro.html">Introduction to <code class="docutils literal notranslate"><span class="pre">quippy</span></code></a></li>
<li class="toctree-l1"><a class="reference internal" href="../install.html">Installation of QUIP and quippy</a></li>
<li class="toctree-l1 current"><a class="reference internal" href="index.html">Tutorials</a><ul class="current">
<li class="toctree-l2"><a class="reference internal" href="Introduction.html">Getting started with <code class="docutils literal notranslate"><span class="pre">Atoms</span></code></a></li>
<li class="toctree-l2"><a class="reference internal" href="quippy-ase-interoperability.html">Interoperability with Atomic Simulation Environment</a></li>
<li class="toctree-l2"><a class="reference internal" href="tutorial.html">Introductory Tutorial</a></li>
<li class="toctree-l2"><a class="reference internal" href="quippy-descriptor-tutorial.html">Using Quippy to compute descriptors of atomic environments</a></li>
<li class="toctree-l2"><a class="reference internal" href="adglass.html">Molecular Dynamics Simulation of Fracture in Quartz</a></li>
<li class="toctree-l2 current"><a class="reference internal" href="adaptive-qmmm.html">Adaptive QM/MM MD of Fracture in Silicon</a><ul class="current">
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-step0.html">Introduction</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-theory.html">Theoretical background</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-step1.html">Step 1: Setup of the Silicon model system</a></li>
<li class="toctree-l3 current"><a class="current reference internal" href="#">Step 2: Classical MD simulation of fracture in Si</a><ul>
<li class="toctree-l4"><a class="reference internal" href="#initialisation-of-the-atomic-system-20-minutes">2.1 Initialisation of the atomic system (20 minutes)</a></li>
<li class="toctree-l4"><a class="reference internal" href="#setup-and-run-the-classical-md-20-minutes">2.2 Setup and run the classical MD (20 minutes)</a></li>
<li class="toctree-l4"><a class="reference internal" href="#visualisation-and-analysis-as-time-permits">2.3 Visualisation and Analysis (as time permits)</a></li>
</ul>
</li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-step3.html">Step 3: LOTF hybrid MD simulation of fracture in Si</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-solutions.html">Solutions</a></li>
<li class="toctree-l3"><a class="reference internal" href="adaptive-qmmm-references.html">References</a></li>
</ul>
</li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../quippy.html">Quippy library reference</a></li>
<li class="toctree-l1"><a class="reference internal" href="../visualisation.html">Visualisation Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../fortran_wrapper.html">Appendix: wrapping Fortran 90 Code</a></li>
<li class="toctree-l1"><a class="reference internal" href="../teach_sparse.html">Command line options of the teach_sparse main program</a></li>
</ul>

            
          
    <p></p>
    <ul>
        <li><a href="../genindex.html">Index</a></li>
        <li><a href="../py-modindex.html">Module Index</a></li>
        <li><a href="../search.html">Search Page</a></li>
        <li><a href="../_modules/index.html">Module Source Code</a></li>
    </ul>
  
        </div>
      </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">

      
      <nav class="wy-nav-top" aria-label="top navigation">
        
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="../index.html">quippy</a>
        
      </nav>


      <div class="wy-nav-content">
        
        <div class="rst-content">
        
          















<div role="navigation" aria-label="breadcrumbs navigation">

  <ul class="wy-breadcrumbs">
    
      <li><a href="../index.html">Docs</a> &raquo;</li>
        
          <li><a href="index.html">Tutorials</a> &raquo;</li>
        
          <li><a href="adaptive-qmmm.html">Adaptive QM/MM MD of Fracture in Silicon</a> &raquo;</li>
        
      <li>Step 2: Classical MD simulation of fracture in Si</li>
    
    
      <li class="wy-breadcrumbs-aside">
        
            
            <a href="../_sources/Tutorials/adaptive-qmmm-step2.rst.txt" rel="nofollow"> View page source</a>
          
        
      </li>
    
  </ul>

  
  <hr/>
</div>
          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
           <div itemprop="articleBody">
            
  
<style>
/* CSS overrides for sphinx_rtd_theme */

/* 24px margin */
.nbinput.nblast,
.nboutput.nblast {
    margin-bottom: 19px;  /* padding has already 5px */
}

/* ... except between code cells! */
.nblast + .nbinput {
    margin-top: -19px;
}

.admonition > p:before {
    margin-right: 4px;  /* make room for the exclamation icon */
}

/* Fix math alignment, see https://github.com/rtfd/sphinx_rtd_theme/pull/686 */
.math {
    text-align: unset;
}
</style>
<div class="section" id="step-2-classical-md-simulation-of-fracture-in-si">
<span id="step2"></span><h1>Step 2: Classical MD simulation of fracture in Si<a class="headerlink" href="#step-2-classical-md-simulation-of-fracture-in-si" title="Permalink to this headline">¶</a></h1>
<p>This part of the tutorial requires the initial crack structure <code class="docutils literal notranslate"><span class="pre">crack.xyz</span></code>
produced in <a class="reference internal" href="adaptive-qmmm-step1.html#step1"><span class="std std-ref">Step 1: Setup of the Silicon model system</span></a>. If you had problems completing the first part, you
can <code class="xref download docutils literal notranslate"><span class="pre">download</span> <span class="pre">it</span> <span class="pre">here</span></code> instead. Start by creating a new
empty script, named <code class="docutils literal notranslate"><span class="pre">run_crack_classical.py</span></code>.</p>
<div class="section" id="initialisation-of-the-atomic-system-20-minutes">
<h2>2.1 Initialisation of the atomic system (20 minutes)<a class="headerlink" href="#initialisation-of-the-atomic-system-20-minutes" title="Permalink to this headline">¶</a></h2>
<div class="section" id="import-the-relevant-modules">
<h3>Import the relevant modules<a class="headerlink" href="#import-the-relevant-modules" title="Permalink to this headline">¶</a></h3>
<p>As before, we import the necessary modules, classes and functions:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>

<span class="kn">from</span> <span class="nn">ase.constraints</span> <span class="k">import</span> <span class="n">FixAtoms</span>
<span class="kn">from</span> <span class="nn">ase.md.verlet</span> <span class="k">import</span> <span class="n">VelocityVerlet</span>
<span class="kn">from</span> <span class="nn">ase.md.velocitydistribution</span> <span class="k">import</span> <span class="n">MaxwellBoltzmannDistribution</span>
<span class="kn">import</span> <span class="nn">ase.units</span> <span class="k">as</span> <span class="nn">units</span>

<span class="kn">from</span> <span class="nn">quippy.atoms</span> <span class="k">import</span> <span class="n">Atoms</span>
<span class="kn">from</span> <span class="nn">quippy.potential</span> <span class="k">import</span> <span class="n">Potential</span>
<span class="kn">from</span> <span class="nn">quippy.io</span> <span class="k">import</span> <span class="n">AtomsWriter</span>

<span class="kn">from</span> <span class="nn">quippy.crack</span> <span class="k">import</span> <span class="p">(</span><span class="n">get_strain</span><span class="p">,</span>
                          <span class="n">get_energy_release_rate</span><span class="p">,</span>
                          <span class="n">ConstantStrainRate</span><span class="p">,</span>
                          <span class="n">find_crack_tip_stress_field</span><span class="p">)</span>
</pre></div>
</div>
<p>Note that the molecular dynamics classes and functions come from
<cite>ASE</cite>, while the <a class="reference internal" href="../atoms.html#quippy.atoms.Atoms" title="quippy.atoms.Atoms"><code class="xref py py-class docutils literal notranslate"><span class="pre">Atoms</span></code></a> object, potentials and
crack-specific functionality come from <cite>quippy</cite>.</p>
</div>
<div class="section" id="definition-of-the-simulation-parameters">
<h3>Definition of the simulation parameters<a class="headerlink" href="#definition-of-the-simulation-parameters" title="Permalink to this headline">¶</a></h3>
<p id="parameters2">Let’s define the simulation parameters. The meaning of each parameter is explained
in the comments on the right of each line:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">input_file</span> <span class="o">=</span> <span class="s1">&#39;crack.xyz&#39;</span>         <span class="c1"># File from which to read crack slab structure</span>
<span class="n">sim_T</span> <span class="o">=</span> <span class="mf">300.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">kB</span>           <span class="c1"># Simulation temperature</span>
<span class="n">nsteps</span> <span class="o">=</span> <span class="mi">10000</span>                   <span class="c1"># Total number of timesteps to run for</span>
<span class="n">timestep</span> <span class="o">=</span> <span class="mf">1.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span>          <span class="c1"># Timestep (NB: time base units are not fs!)</span>
<span class="n">cutoff_skin</span> <span class="o">=</span> <span class="mf">2.0</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">Ang</span>      <span class="c1"># Amount by which potential cutoff is increased</span>
                                 <span class="c1"># for neighbour calculations</span>
<span class="n">tip_move_tol</span> <span class="o">=</span> <span class="mf">10.0</span>              <span class="c1"># Distance tip has to move before crack</span>
                                 <span class="c1"># is taken to be running</span>
<span class="n">strain_rate</span> <span class="o">=</span> <span class="mf">1e-5</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">fs</span><span class="p">)</span>  <span class="c1"># Strain rate</span>
<span class="n">traj_file</span> <span class="o">=</span> <span class="s1">&#39;traj.nc&#39;</span>            <span class="c1"># Trajectory output file (NetCDF format)</span>
<span class="n">print_interval</span> <span class="o">=</span> <span class="mi">10</span>              <span class="c1"># Number of time steps between</span>
                                 <span class="c1"># writing output frames</span>
<span class="n">param_file</span> <span class="o">=</span> <span class="s1">&#39;params.xml&#39;</span>        <span class="c1"># Filename of XML file containing</span>
                                 <span class="c1"># potential parameters</span>
<span class="n">mm_init_args</span> <span class="o">=</span> <span class="s1">&#39;IP SW&#39;</span>           <span class="c1"># Initialisation arguments for</span>
                                 <span class="c1"># classical potential</span>
</pre></div>
</div>
</div>
<div class="section" id="setup-of-the-atomic-structure">
<h3>Setup of the atomic structure<a class="headerlink" href="#setup-of-the-atomic-structure" title="Permalink to this headline">¶</a></h3>
<p>As a first step, we need to initialise the
<a class="reference internal" href="../atoms.html#quippy.atoms.Atoms" title="quippy.atoms.Atoms"><code class="xref py py-class docutils literal notranslate"><span class="pre">Atoms</span></code></a> object by loading the atomic structure created
in <a class="reference internal" href="adaptive-qmmm-step1.html#step1"><span class="std std-ref">Step 1: Setup of the Silicon model system</span></a> from the <cite>input_file</cite> <code class="docutils literal notranslate"><span class="pre">crack.xyz</span></code>. Note that the global
<cite>fortran_indexing</cite> setting should be set to <cite>False</cite>. Otherwise <cite>quippy</cite> uses
atom indices in the range <span class="math notranslate nohighlight">\(1 \ldots N\)</span>, which would not be consistent with
the python indexing used in ASE (<span class="math notranslate nohighlight">\(0\ldots N-1\)</span>).</p>
<p>It is also necessary to read in the original width and height of the slab and
the original crack position, which were saved in <code class="docutils literal notranslate"><span class="pre">crack.xyz</span></code> at the end
of <a class="reference internal" href="adaptive-qmmm-step1.html#step1"><span class="std std-ref">Step 1</span></a>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="nb">print</span> <span class="s1">&#39;Loading atoms from file </span><span class="si">%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">input_file</span>
<span class="n">atoms</span> <span class="o">=</span> <span class="o">...</span>                                     <span class="c1"># Load atoms from file</span>

<span class="n">orig_height</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span>          <span class="c1"># Initialise original height</span>
<span class="n">orig_crack_pos</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">]</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>  <span class="c1"># Initialise original crack position</span>
</pre></div>
</div>
<p>Note that we make a copy of the <code class="docutils literal notranslate"><span class="pre">CrackPos</span></code> entry in the
<code class="xref py py-attr docutils literal notranslate"><span class="pre">info</span></code> dictionary, since otherwise
<cite>orig_crack_pos</cite> would continue to refer to the current crack position
as it is updated during the dynamical simulation.</p>
</div>
<div class="section" id="setup-of-the-constraints">
<h3>Setup of the constraints<a class="headerlink" href="#setup-of-the-constraints" title="Permalink to this headline">¶</a></h3>
<p>Now we can set constraints on the atomic structure which will be
enforced during dynamics. In order to carry out the fracture MD
simulation, we need to fix the positions of the top and bottom atomic
rows (we call this constraint <cite>fix_atoms</cite>).</p>
<p>The <cite>fix_atoms</cite> constraint, which is exactly the same as the constraint used
for <a class="reference internal" href="adaptive-qmmm-step1.html#crack-fixatoms"><span class="std std-ref">relaxing the positions of the crack slab</span></a> above. In
order to do this, we need to find the <cite>y</cite> coordinate of the top, bottom atomic
rows. The <cite>x</cite> coordinates of the left and right edges of the slab might also
be useful later on. This can be easily done as before:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">top</span> <span class="o">=</span> <span class="o">...</span>     <span class="c1"># Maximum y coordinate</span>
<span class="n">bottom</span> <span class="o">=</span> <span class="o">...</span>  <span class="c1"># Minimum y coordinate</span>
<span class="n">left</span> <span class="o">=</span> <span class="o">...</span>    <span class="c1"># Minimum x coordinate</span>
<span class="n">right</span> <span class="o">=</span> <span class="o">...</span>   <span class="c1"># Maximum x coordinate</span>
</pre></div>
</div>
<p>Now it is possible to define the <cite>fixed_mask</cite> array, which is <cite>True</cite>
for each atom whose position needs to be fixed, and <cite>False</cite> otherwise,
exactly as before, and to initialise the <cite>fix_atoms</cite> constraint, in
the same way we did it in <cite>Step 1</cite> (i.e., using the
<code class="xref py py-class docutils literal notranslate"><span class="pre">FixAtoms</span></code> class):</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">fixed_mask</span> <span class="o">=</span> <span class="o">...</span>                             <span class="c1"># Define the list of fixed atoms</span>
<span class="n">fix_atoms</span> <span class="o">=</span> <span class="o">...</span>                              <span class="c1"># Initialise the constraint</span>
<span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Fixed </span><span class="si">%d</span><span class="s1"> atoms</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">fixed_mask</span><span class="o">.</span><span class="n">sum</span><span class="p">())</span> <span class="c1"># Print the number of fixed atoms</span>
</pre></div>
</div>
<p>This constraint needs to be attached to our <cite>atoms</cite> object
(see <a class="reference internal" href="../atoms.html#quippy.atoms.Atoms.set_constraint" title="quippy.atoms.Atoms.set_constraint"><code class="xref py py-meth docutils literal notranslate"><span class="pre">set_constraint()</span></code></a>):</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">atoms</span><span class="o">.</span> <span class="o">...</span>  <span class="c1"># Attach the constraints to atoms</span>
</pre></div>
</div>
<p>To increase <span class="math notranslate nohighlight">\(\epsilon_{yy}\)</span> of all atoms at a constant rate (see
the <cite>strain_rate</cite> and <cite>timestep</cite> <a class="reference internal" href="#parameters2"><span class="std std-ref">parameters</span></a>), we
use the <a class="reference internal" href="../crack.html#quippy.crack.ConstantStrainRate" title="quippy.crack.ConstantStrainRate"><code class="xref py py-class docutils literal notranslate"><span class="pre">ConstantStrainRate</span></code></a> class:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">strain_atoms</span> <span class="o">=</span> <span class="n">ConstantStrainRate</span><span class="p">(</span><span class="n">orig_height</span><span class="p">,</span> <span class="n">strain_rate</span><span class="o">*</span><span class="n">timestep</span><span class="p">)</span>
</pre></div>
</div>
<p>You can look at the documentation for the
<a class="reference internal" href="../crack.html#quippy.crack.ConstantStrainRate" title="quippy.crack.ConstantStrainRate"><code class="xref py py-class docutils literal notranslate"><span class="pre">ConstantStrainRate</span></code></a> class to see how this works. The
<code class="xref py py-meth docutils literal notranslate"><span class="pre">adjust_positions()</span></code> routine simply increases the strain of
all atoms. We will attach this to the dynamics in step 2.2 below.</p>
</div>
<div class="section" id="setup-of-the-potential">
<h3>Setup of the potential<a class="headerlink" href="#setup-of-the-potential" title="Permalink to this headline">¶</a></h3>
<p>Before starting the MD simulation, the SW classical potential must be
initialised and attached to the <cite>atoms</cite> object.  As in <cite>Step 1</cite>, we
use quippy’s <a class="reference internal" href="../potential.html#quippy.potential.Potential" title="quippy.potential.Potential"><code class="xref py py-class docutils literal notranslate"><span class="pre">Potential</span></code></a> class, but now we
need to pass the <cite>cutoff_skin</cite> parameter, which is used to decide when
the neighbour list needs to be updated (see the attribute
<code class="xref py py-attr docutils literal notranslate"><span class="pre">cutoff_skin</span></code>). Moreover, we request
the potential to compute per-atom stresses whenever we compute forces
using <code class="xref py py-meth docutils literal notranslate"><span class="pre">set_default_quantities()</span></code>, to
save time when locating the crack tip (discussed in more detail
<a class="reference internal" href="#position-crack-tip"><span class="std std-ref">below</span></a>). The
<a class="reference internal" href="../atoms.html#quippy.atoms.Atoms.set_calculator" title="quippy.atoms.Atoms.set_calculator"><code class="xref py py-meth docutils literal notranslate"><span class="pre">set_calculator()</span></code></a> method should then be used
to set the calculator to the SW potential:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">mm_pot</span> <span class="o">=</span> <span class="o">...</span>   <span class="c1"># Initialise the SW potential with cutoff_skin</span>
<span class="n">mm_pot</span><span class="o">.</span><span class="n">set_default_quantities</span><span class="p">([</span><span class="s1">&#39;stresses&#39;</span><span class="p">])</span>
<span class="n">atoms</span><span class="o">.</span> <span class="o">...</span>     <span class="c1"># Set the calculator</span>
</pre></div>
</div>
</div>
<div class="section" id="milestone-2-1">
<h3>Milestone 2.1<a class="headerlink" href="#milestone-2-1" title="Permalink to this headline">¶</a></h3>
<p>At this stage your script should look something like <a class="reference download internal" download="" href="../_downloads/88e6e0210271f9314e0fe935df418a02/run_crack_classical_1.py"><code class="xref download docutils literal notranslate"><span class="pre">this</span></code></a>.</p>
</div>
</div>
<div class="section" id="setup-and-run-the-classical-md-20-minutes">
<h2>2.2 Setup and run the classical MD (20 minutes)<a class="headerlink" href="#setup-and-run-the-classical-md-20-minutes" title="Permalink to this headline">¶</a></h2>
<div class="section" id="setting-initial-velocities-and-constructing-the-dynamics-object">
<h3>Setting initial velocities and constructing the dynamics object<a class="headerlink" href="#setting-initial-velocities-and-constructing-the-dynamics-object" title="Permalink to this headline">¶</a></h3>
<p>There are still a few things that need to be done before running the
MD fracture simulation. We will follow the standard <a class="reference external" href="https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html">ASE molecular
dynamics methodology</a>. We will set
the initial temperature of the system to <cite>2*sim_T</cite>: it will then
equilibrate to <cite>sim_T</cite>, by the Virial theorem:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">MaxwellBoltzmannDistribution</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="mf">2.0</span><span class="o">*</span><span class="n">sim_T</span><span class="p">)</span>
</pre></div>
</div>
<p>A MD simulation in the NVE ensemble, using the Velocity Verlet
algorithm, can be initialised with the ASE
<code class="xref py py-class docutils literal notranslate"><span class="pre">VelocityVerlet</span></code> class, which requires two
arguments: the atoms and the time step (which should come from the
<cite>timestep</cite> <a class="reference internal" href="#parameters2"><span class="std std-ref">parameter</span></a>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">dynamics</span> <span class="o">=</span> <span class="o">...</span>   <span class="c1"># Initialise the dynamics</span>
</pre></div>
</div>
</div>
<div class="section" id="printing-status-information">
<h3>Printing status information<a class="headerlink" href="#printing-status-information" title="Permalink to this headline">¶</a></h3>
<p>Let’s also define a function that prints the relevant information at
each time step of the MD simulation. The information can be saved
inside the <a class="reference internal" href="../atoms.html#quippy.atoms.Atoms.info" title="quippy.atoms.Atoms.info"><code class="xref py py-attr docutils literal notranslate"><span class="pre">info</span></code></a> dictionary, so that it
also gets saved to the trajectory file <cite>traj_file</cite>.</p>
<p>The elapsed simulation time can be obtained with
<code class="docutils literal notranslate"><span class="pre">dynamics.get_time()</span></code> (note that the time unit in ASE is
<span class="math notranslate nohighlight">\(\mathrm{\AA}\sqrt{\mathrm{amu}/\mathrm{eV}}\)</span>, not <cite>fs</cite>). You
should use the <code class="xref py py-meth docutils literal notranslate"><span class="pre">get_kinetic_energy()</span></code> method to
calculate the temperature (<em>Note</em>: you will need the <code class="xref py py-attr docutils literal notranslate"><span class="pre">units.kB</span></code>
constant, which gives the value of the Boltzmann constant in eV/K),
and the functions <a class="reference internal" href="../crack.html#quippy.crack.get_strain" title="quippy.crack.get_strain"><code class="xref py py-func docutils literal notranslate"><span class="pre">get_strain()</span></code></a> and
<a class="reference internal" href="../crack.html#quippy.crack.get_energy_release_rate" title="quippy.crack.get_energy_release_rate"><code class="xref py py-func docutils literal notranslate"><span class="pre">get_energy_release_rate()</span></code></a> to return the current
strain energy release rate, respectively.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">printstatus</span><span class="p">():</span>
    <span class="k">if</span> <span class="n">dynamics</span><span class="o">.</span><span class="n">nsteps</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="nb">print</span> <span class="s2">&quot;&quot;&quot;</span>
<span class="s2">State      Time/fs    Temp/K     Strain      G/(J/m^2)  CrackPos/A D(CrackPos)/A</span>
<span class="s2">---------------------------------------------------------------------------------&quot;&quot;&quot;</span>

    <span class="n">log_format</span> <span class="o">=</span> <span class="p">(</span><span class="s1">&#39;</span><span class="si">%(label)-4s%(time)12.1f%(temperature)12.6f</span><span class="s1">&#39;</span><span class="o">+</span>
        <span class="s1">&#39;</span><span class="si">%(strain)12.5f%(G)12.4f%(crack_pos_x)12.2f</span><span class="s1">    (</span><span class="si">%(d_crack_pos_x)+5.2f</span><span class="s1">)&#39;</span><span class="p">)</span>

    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;label&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="s1">&#39;D&#39;</span>                <span class="c1"># Label for the status line</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;time&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="o">...</span>                 <span class="c1"># Get simulation time</span>
                                             <span class="c1"># and convert to fs</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;temperature&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="o">...</span>          <span class="c1"># Get temperature in K</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;strain&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="o">...</span>               <span class="c1"># Get strain</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="o">...</span>                    <span class="c1"># Get energy release rate,</span>
                                             <span class="c1"># and convert to J/m^2</span>
    <span class="n">crack_pos</span> <span class="o">=</span> <span class="o">...</span>                          <span class="c1"># Find crack tip as in step 1</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;crack_pos_x&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;d_crack_pos_x&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">orig_crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

    <span class="nb">print</span> <span class="n">log_format</span> <span class="o">%</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span>
</pre></div>
</div>
<p>This logger can be now attached to the <cite>dynamics</cite>, so that the information is
printed at every time step during the simulations:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">printstatus</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="checking-if-the-crack-has-advanced">
<h3>Checking if the crack has advanced<a class="headerlink" href="#checking-if-the-crack-has-advanced" title="Permalink to this headline">¶</a></h3>
<p>The same can be done to check during the simulation if the crack has advanced,
and to stop incrementing the strain if it has:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">check_if_cracked</span><span class="p">(</span><span class="n">atoms</span><span class="p">):</span>
    <span class="n">crack_pos</span> <span class="o">=</span> <span class="o">...</span>                          <span class="c1"># Find crack tip position</span>

    <span class="c1"># stop straining if crack has advanced more than tip_move_tol</span>
    <span class="k">if</span> <span class="ow">not</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]</span> <span class="ow">and</span> <span class="p">(</span><span class="n">crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">orig_crack_pos</span><span class="p">[</span><span class="mi">0</span><span class="p">])</span> <span class="o">&gt;</span> <span class="n">tip_move_tol</span><span class="p">:</span>
        <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;is_cracked&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="k">del</span> <span class="n">atoms</span><span class="o">.</span><span class="n">constraints</span><span class="p">[</span><span class="n">atoms</span><span class="o">.</span><span class="n">constraints</span><span class="o">.</span><span class="n">index</span><span class="p">(</span><span class="n">strain_atoms</span><span class="p">)]</span>
</pre></div>
</div>
<p>The <cite>check_if_cracked</cite> function can now be attached to the dynamical
system, requesting an interval of 1 step (i.e. every time) and passing the
<cite>atoms</cite> object along to the function as an extra argument:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">check_if_cracked</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">atoms</span><span class="p">)</span>
</pre></div>
</div>
<p>We also need to attach the <cite>:meth:`quippy.crack.ConstrainStrainRate.apply_strain</cite> method
of <cite>strain_atoms</cite> to the dynamics:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">dynamics</span><span class="o">.</span><span class="n">attach</span><span class="p">(</span><span class="n">strain_atoms</span><span class="o">.</span><span class="n">apply_strain</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="n">atoms</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="saving-the-trajectory">
<h3>Saving the trajectory<a class="headerlink" href="#saving-the-trajectory" title="Permalink to this headline">¶</a></h3>
<p>Finally, we need to initialise the trajectory file <cite>traj_file</cite> and to
save frames to the trajectory every <cite>traj_interval</cite> time steps. This
is done by creating a trajectory object with the
<a class="reference internal" href="../io.html#quippy.io.AtomsWriter" title="quippy.io.AtomsWriter"><code class="xref py py-func docutils literal notranslate"><span class="pre">AtomsWriter()</span></code></a> function, and then attaching this
trajectory to the <cite>dynamics</cite>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">trajectory</span> <span class="o">=</span> <span class="o">...</span>    <span class="c1"># Initialise the trajectory</span>
<span class="n">dynamics</span><span class="o">.</span> <span class="o">...</span>       <span class="c1"># Attach the trajectory with an interval of</span>
                    <span class="c1"># traj_interval, passing atoms as an extra argument</span>
</pre></div>
</div>
<p>We will save the trajectory in <a class="reference internal" href="../io.html#netcdf"><span class="std std-ref">NetCDF</span></a> format. This is a binary
file format that is similar with the <a class="reference internal" href="../io.html#extendedxyz"><span class="std std-ref">Extended XYZ</span></a> format we used
earlier, with the advantage of being more efficient for large files.</p>
</div>
<div class="section" id="running-the-dynamics">
<h3>Running the dynamics<a class="headerlink" href="#running-the-dynamics" title="Permalink to this headline">¶</a></h3>
<p>After all this, a single command will run the MD for <cite>nsteps</cite> (see the <a class="reference external" href="https://wiki.fysik.dtu.dk/ase/tutorials/md/md.html">ASE
molecular dynamics methodology</a> for more information):</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">dynamics</span><span class="o">.</span><span class="n">run</span><span class="p">(</span><span class="n">nsteps</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="section" id="milestone-2-2">
<h3>Milestone 2.2<a class="headerlink" href="#milestone-2-2" title="Permalink to this headline">¶</a></h3>
<p>If you have problems you can look at the complete version of the
<a class="reference internal" href="adaptive-qmmm-solutions.html#run-crack-classical"><span class="std std-ref">Step 2 solution — run_crack_classical.py</span></a> script. Leave your classical MD simulation
running and move on to the next section of the tutorial.</p>
<p>The first few lines produced by the <code class="docutils literal notranslate"><span class="pre">run_crack_classical.py</span></code> script should
look something like this:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">Loading</span> <span class="n">atoms</span> <span class="kn">from</span> <span class="nn">file</span> <span class="n">crack</span><span class="o">.</span><span class="n">xyz</span>
<span class="n">Fixed</span> <span class="mi">240</span> <span class="n">atoms</span>


<span class="n">State</span>      <span class="n">Time</span><span class="o">/</span><span class="n">fs</span>    <span class="n">Temp</span><span class="o">/</span><span class="n">K</span>     <span class="n">Strain</span>      <span class="n">G</span><span class="o">/</span><span class="p">(</span><span class="n">J</span><span class="o">/</span><span class="n">m</span><span class="o">^</span><span class="mi">2</span><span class="p">)</span>  <span class="n">CrackPos</span><span class="o">/</span><span class="n">A</span> <span class="n">D</span><span class="p">(</span><span class="n">CrackPos</span><span class="p">)</span><span class="o">/</span><span class="n">A</span>
<span class="o">---------------------------------------------------------------------------------</span>
<span class="n">D</span>            <span class="mf">1.0</span>  <span class="mf">560.097755</span>     <span class="mf">0.08427</span>      <span class="mf">5.0012</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.00</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">2.0</span>  <span class="mf">550.752265</span>     <span class="mf">0.08428</span>      <span class="mf">5.0024</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.00</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">3.0</span>  <span class="mf">535.568949</span>     <span class="mf">0.08429</span>      <span class="mf">5.0036</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.00</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">4.0</span>  <span class="mf">515.074874</span>     <span class="mf">0.08430</span>      <span class="mf">5.0047</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.00</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">5.0</span>  <span class="mf">489.977973</span>     <span class="mf">0.08431</span>      <span class="mf">5.0059</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.00</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">6.0</span>  <span class="mf">461.140488</span>     <span class="mf">0.08432</span>      <span class="mf">5.0071</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.00</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">7.0</span>  <span class="mf">429.546498</span>     <span class="mf">0.08433</span>      <span class="mf">5.0083</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.00</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">8.0</span>  <span class="mf">396.264666</span>     <span class="mf">0.08434</span>      <span class="mf">5.0095</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.01</span><span class="p">)</span>
<span class="n">D</span>            <span class="mf">9.0</span>  <span class="mf">362.407525</span>     <span class="mf">0.08435</span>      <span class="mf">5.0107</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.01</span><span class="p">)</span>
<span class="n">D</span>           <span class="mf">10.0</span>  <span class="mf">329.088872</span>     <span class="mf">0.08436</span>      <span class="mf">5.0119</span>      <span class="o">-</span><span class="mf">30.61</span>    <span class="p">(</span><span class="o">-</span><span class="mf">0.01</span><span class="p">)</span>
</pre></div>
</div>
<p>Here we see the current time, temperature, strain, energy release rate <cite>G</cite>, the
<cite>x</cite> coordinate of the crack position, and the change in the crack position since
the beginning of the simulation. In the early stages of the calculation, the
strain and <cite>G</cite> are both increasing, and the temperature is rapidly falling
towards <code class="docutils literal notranslate"><span class="pre">sim_T</span> <span class="pre">=</span> <span class="pre">300</span></code> as anticipated.</p>
</div>
</div>
<div class="section" id="visualisation-and-analysis-as-time-permits">
<span id="visualisation2"></span><h2>2.3 Visualisation and Analysis (as time permits)<a class="headerlink" href="#visualisation-and-analysis-as-time-permits" title="Permalink to this headline">¶</a></h2>
<p>Start another <cite>ipython</cite> session is a new terminal with plotting
support enabled, using the shell command:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">ipython</span> <span class="o">--</span><span class="n">pylab</span>
</pre></div>
</div>
<p>This will allow you to look at the progress of your classical fracture
simulation while it continues to run. All the example code given in
this section should be entered directly at the <cite>ipython</cite> prompt.</p>
<p>The first step is to import everything from <cite>quippy</cite> using the
<a class="reference internal" href="../qlab.html#module-qlab" title="qlab: Interactive visualisation of atomic systems"><code class="xref py py-mod docutils literal notranslate"><span class="pre">qlab</span></code></a> interactive module, then open your trajectory using the
<a class="reference internal" href="../qlab.html#qlab.view" title="qlab.view"><code class="xref py py-func docutils literal notranslate"><span class="pre">view()</span></code></a> function:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="kn">from</span> <span class="nn">qlab</span> <span class="k">import</span> <span class="o">*</span>
<span class="n">set_fortran_indexing</span><span class="p">(</span><span class="kc">False</span><span class="p">)</span>
<span class="n">view</span><span class="p">(</span><span class="s2">&quot;traj.nc&quot;</span><span class="p">)</span>
</pre></div>
</div>
<p>As we saw <a class="reference internal" href="adaptive-qmmm-step1.html#latticeconstant"><span class="std std-ref">earlier</span></a>, this will open an AtomEye viewer
window containing a visual representation of your crack system (as before
<code class="docutils literal notranslate"><span class="pre">fortran_indexing=False</span></code> is used to number the atoms starting from zero). You
can use the <cite>Insert</cite> and <cite>Delete</cite> keys to move forwards or backwards through the
trajectory, or <cite>Ctrl+Insert</cite> and <cite>Ctrl+Delete</cite> to jump to the first or last
frame – note that the focus must be on the AtomEye viewer window when you use
any keyboard shortcuts. The current frame number is shown in the title bar of
the window.</p>
<p>The function <a class="reference internal" href="../qlab.html#qlab.gcat" title="qlab.gcat"><code class="xref py py-func docutils literal notranslate"><span class="pre">gcat()</span></code></a> (short for “get current atoms”) returns a
reference to the <a class="reference internal" href="../atoms.html#quippy.atoms.Atoms" title="quippy.atoms.Atoms"><code class="xref py py-class docutils literal notranslate"><span class="pre">Atoms</span></code></a> object currently being visualised (i.e. to the
current frame from the trajectory file). Similarly, the <a class="reference internal" href="../qlab.html#qlab.gcv" title="qlab.gcv"><code class="xref py py-func docutils literal notranslate"><span class="pre">gcv()</span></code></a>
function returns a reference to the entire trajectory currently being viewed as
an <a class="reference internal" href="../qlab.html#qlab.AtomsReaderViewer" title="qlab.AtomsReaderViewer"><code class="xref py py-class docutils literal notranslate"><span class="pre">AtomsReaderViewer</span></code></a> object.</p>
<p>You can change the frame increment rate by setting
the <a class="reference internal" href="../atomeye.html#atomeye.AtomEyeViewer.delta" title="atomeye.AtomEyeViewer.delta"><code class="xref py py-attr docutils literal notranslate"><span class="pre">delta</span></code></a> attribute of the viewer, e.g. to
advance by ten frames at a time:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">set_delta</span><span class="p">(</span><span class="mi">10</span><span class="p">)</span>
</pre></div>
</div>
<p>Or, to jump directly to frame 100:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">set_frame</span><span class="p">(</span><span class="mi">100</span><span class="p">)</span>
</pre></div>
</div>
<p>You can repeat the <code class="docutils literal notranslate"><span class="pre">view(&quot;traj.nc&quot;)</span></code>
command as your simulation progresses to reload the file (you can use <cite>Ctrl+R</cite>
in the <cite>ipython</cite> console to search backwards in the session history to save
typing).</p>
<div class="section" id="stress-field-analysis">
<span id="stress-analysis"></span><h3>Stress field analysis<a class="headerlink" href="#stress-field-analysis" title="Permalink to this headline">¶</a></h3>
<p>To compute and display the instantaneous principal per-atom stress
<span class="math notranslate nohighlight">\(\sigma_{yy}\)</span> as computed by the SW potential for a
configuration near the beginning of your dynamical simulation:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">mm_pot</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="s1">&#39;IP SW&#39;</span><span class="p">,</span> <span class="n">param_filename</span><span class="o">=</span><span class="s1">&#39;params.xml&#39;</span><span class="p">)</span>
<span class="n">at</span> <span class="o">=</span> <span class="n">gcat</span><span class="p">()</span>
<span class="n">at</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">mm_pot</span><span class="p">)</span>
<span class="n">mm_sigma</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">get_stresses</span><span class="p">()</span>
<span class="n">sigma_yy</span> <span class="o">=</span> <span class="n">mm_sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span>
<span class="n">aux_property_coloring</span><span class="p">(</span><span class="n">sigma_yy</span><span class="p">)</span>
</pre></div>
</div>
<p>The <cite>mm_sigma</cite> array has shape <cite>(len(atoms), 3, 3)</cite>, i.e. it is
made up of a <span class="math notranslate nohighlight">\(3 \times 3\)</span> stress tensor <span class="math notranslate nohighlight">\(\sigma_{ij}\)</span> for
each atom. The <cite>sigma_yy</cite> array is the <code class="docutils literal notranslate"><span class="pre">[1,</span> <span class="pre">1]</span></code> component of each of
these arrays, i.e. <span class="math notranslate nohighlight">\(\sigma_{yy}\)</span>. To read off the value of the
stress on a particular atom, just <cite>right click</cite> on it. As before, this
prints various information in the <cite>ipython</cite> console. The <cite>_show</cite>
property corresponds to the values currently being used to colour the
atoms. You will see that <span class="math notranslate nohighlight">\(\sigma_{yy}\)</span> is very strongly peaked
near the crack tip. If you prefer to see the values in GPa, you could
do</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">aux_property_coloring</span><span class="p">(</span><span class="n">sigma_yy</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">GPa</span><span class="p">)</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/sigma_yy.png"><img alt="../_images/sigma_yy.png" class="align-center" src="../_images/sigma_yy.png" style="width: 600px;" /></a>
<p>The concept of per-atom stresses is a little arbitrary. The values we
are plotting here were obtained from partitioning the total virial
stress tensor, which is given by</p>
<div class="math notranslate nohighlight">
\[\begin{split}\tau_{ij} = \frac{1}{\Omega} \sum_{k \in \Omega} (-m^{(k)} (u_i^{(k)}-
\bar{u}_i) (u_j^{(k)}- \bar{u}_j) %\\
 +  \frac{1}{2} \sum_{\ell \in \Omega} ( x_i^{(\ell)} - x_i^{(k)}) f_j^{(k\ell)}
)\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(k\)</span> and <span class="math notranslate nohighlight">\(l\)</span> are atom indices, <span class="math notranslate nohighlight">\(ijk\)</span> are Cartesian
indicies, <span class="math notranslate nohighlight">\(\Omega\)</span> is the cell volume, <span class="math notranslate nohighlight">\(m^{(k)}\)</span>,
<span class="math notranslate nohighlight">\(u^{(k)}\)</span>, <span class="math notranslate nohighlight">\(x^{(k)}\)</span> and <span class="math notranslate nohighlight">\(f^{(k)}\)</span> are respectively the
mass, velocity, position of atom <span class="math notranslate nohighlight">\(k\)</span> and <span class="math notranslate nohighlight">\(f^{kl}_j\)</span> is
the <span class="math notranslate nohighlight">\(j\)</span>th component of the force between atoms <span class="math notranslate nohighlight">\(k\)</span> and
<span class="math notranslate nohighlight">\(l\)</span>. The first term is a kinetic contribution which vanishes at
near zero temperature, and it is common to use the second term to
define a per-atom stress tensor.</p>
<p>Note, however, that this requires a definition of the atomic volume. By default
the <a class="reference internal" href="../potential.html#quippy.potential.Potential.get_stresses" title="quippy.potential.Potential.get_stresses"><code class="xref py py-meth docutils literal notranslate"><span class="pre">get_stresses()</span></code></a> function simply divides the
total cell volume by the number of atoms to get the volume per atom. This is
not a very good approximation for our cell, which contains a lot of empty
vacuum, so the volume per atom comes out much too large, and the stress
components much too small, e.g. the peak stress, which you can print in units of
GPa with:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="nb">print</span> <span class="n">mm_sigma</span><span class="o">.</span><span class="n">max</span><span class="p">()</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">GPa</span>
</pre></div>
</div>
<p>is around 4 GPa. Values of stress in better agreement with linear
elastic theory can be obtained by assuming all atoms occupy the same
volume as they would in the equilibrium bulk structure:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">mm_pot</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="n">vol_per_atom</span><span class="o">=</span><span class="n">si_bulk</span><span class="o">.</span><span class="n">get_volume</span><span class="p">()</span><span class="o">/</span><span class="nb">len</span><span class="p">(</span><span class="n">si_bulk</span><span class="p">))</span>
<span class="n">mm_sigma</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">get_stresses</span><span class="p">()</span>
<span class="nb">print</span> <span class="n">mm_sigma</span><span class="o">.</span><span class="n">max</span><span class="p">()</span><span class="o">/</span><span class="n">units</span><span class="o">.</span><span class="n">GPa</span>
</pre></div>
</div>
<p>gives a value of around 25 GPa. As this is only a simple rescaling,
the unscaled virial stress values are perfectly adequate for locating
the crack tip.</p>
<p>Use values from the <cite>sigma_yy</cite> array to plot the <span class="math notranslate nohighlight">\(\sigma_{yy}\)</span> virial
stress along the line <span class="math notranslate nohighlight">\(y=0\)</span> ahead of the crack tip, and verify the stress
obeys the expected <span class="math notranslate nohighlight">\(1/\sqrt{r}\)</span> divergence near the crack tip, and tends
to a constant value ahead of the crack, due to the thin strip loading. <em>Hint:</em>
use a mask to select the relevant atoms, as we did when fixing the edge atoms
above. You can use the matplotlib <a class="reference external" href="https://matplotlib.org/api/_as_gen/matplotlib.pyplot.plot.html#matplotlib.pyplot.plot" title="(in Matplotlib v3.0.3)"><code class="xref py py-func docutils literal notranslate"><span class="pre">plot()</span></code></a> function to
produce a plot.</p>
</div>
<div class="section" id="time-averaged-stress-field">
<span id="time-avg-stress"></span><h3>Time-averaged stress field<a class="headerlink" href="#time-averaged-stress-field" title="Permalink to this headline">¶</a></h3>
<p>By now, you should have a few picoseconds of dynamics in your trajectory file.
Reload with <code class="docutils literal notranslate"><span class="pre">view(&quot;traj.nc&quot;)</span></code> to see what is happening. You can jump to the
end with <cite>Ctrl+Delete</cite>, or by typing <cite>last()</cite> into the <cite>ipython</cite> console. Here
is what the instantaneous <span class="math notranslate nohighlight">\(\sigma_{yy}\)</span> looks like after 5 ps of dynamics:</p>
<a class="reference internal image-reference" href="../_images/classical-crack-sigma-yy.png"><img alt="../_images/classical-crack-sigma-yy.png" class="align-center" src="../_images/classical-crack-sigma-yy.png" style="width: 600px;" /></a>
<p>As you can see, the stress field is rather noisy because of
contributions made by the random thermal motion of atoms. The
<a class="reference internal" href="../crack.html#quippy.crack.find_crack_tip_stress_field" title="quippy.crack.find_crack_tip_stress_field"><code class="xref py py-func docutils literal notranslate"><span class="pre">find_crack_tip_stress_field()</span></code></a> uses an exponential
moving average of the stress field when finding the tip. This average
is stored in the <code class="docutils literal notranslate"><span class="pre">avg_sigma</span></code> <code class="xref py py-attr docutils literal notranslate"><span class="pre">array</span> <span class="pre">entry</span></code> inside the Atoms object, which is saved
with each frame in the trajectory. For techical reasons this is stored
as a reshaped array of shape <code class="docutils literal notranslate"><span class="pre">(len(atoms),</span> <span class="pre">9)</span></code> rather than
<code class="docutils literal notranslate"><span class="pre">(len(atoms),</span> <span class="pre">3,</span> <span class="pre">3)</span></code> array, so you can find the <span class="math notranslate nohighlight">\(sigma_{yy}\)</span>
components in the 5th column (counting from zero as usual in Python),
i.e.</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">aux_property_coloring</span><span class="p">(</span><span class="n">gcat</span><span class="p">()</span><span class="o">.</span><span class="n">arrays</span><span class="p">[</span><span class="s1">&#39;avg_sigma&#39;</span><span class="p">][:,</span> <span class="mi">4</span><span class="p">])</span>
</pre></div>
</div>
<p>You should find that the crack tip is more well defined in the average stress:</p>
<a class="reference internal image-reference" href="../_images/classical-crack-sigma-yy-average.png"><img alt="../_images/classical-crack-sigma-yy-average.png" class="align-center" src="../_images/classical-crack-sigma-yy-average.png" style="width: 600px;" /></a>
</div>
<div class="section" id="geometry-and-coordination-analysis">
<span id="coordination"></span><h3>Geometry and coordination analysis<a class="headerlink" href="#geometry-and-coordination-analysis" title="Permalink to this headline">¶</a></h3>
<p>Press <cite>k</cite> to colour the atoms by coordination. This is based on the
<a class="reference internal" href="../atoms.html#quippy.atoms.Atoms.nneightol" title="quippy.atoms.Atoms.nneightol"><code class="xref py py-attr docutils literal notranslate"><span class="pre">nneightol</span></code></a> attribute of the Atoms object, which we set
to a value of 1.3 in the <code class="docutils literal notranslate"><span class="pre">make_crack.py</span></code> script. This factor acts as
a multipler for the covalent radii of the atomic species, taken from
the <code class="xref py py-attr docutils literal notranslate"><span class="pre">quippy.periodictable.ElementCovRad</span></code> array. You can check
the maximum Si–Si bond-length this corresponds to with:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="nb">print</span> <span class="mf">1.3</span><span class="o">*</span><span class="mi">2</span><span class="o">*</span><span class="n">ElementCovRad</span><span class="p">[</span><span class="mi">14</span><span class="p">]</span>
</pre></div>
</div>
<p>Note that <code class="docutils literal notranslate"><span class="pre">14</span></code> is the atomic number of silicon. After the simulation has run
for a little while, you should be able to see both under-coordinated (green) and
over-coordinated (red) atoms near the crack tip.</p>
<p>Here is a typical snapshot at the end of 10 ps of dynamics. Note the
large number of defects, indicating that the fracture surface is not
atomically smooth as we find it to be in experiments. In your
simulation you may be able to spot signs of energy dissipation
mechanisms, such as dislocation emission from the crack tip.</p>
<a class="reference internal image-reference" href="../_images/classical-crack-coordination.png"><img alt="../_images/classical-crack-coordination.png" class="align-center" src="../_images/classical-crack-coordination.png" style="width: 600px;" /></a>
</div>
<div class="section" id="rendering-a-movie-of-the-simulation">
<span id="render-movie"></span><h3>Rendering a movie of the simulation<a class="headerlink" href="#rendering-a-movie-of-the-simulation" title="Permalink to this headline">¶</a></h3>
<p>If you would like to make a movie of your simulation, you can use
the <a class="reference internal" href="../qlab.html#qlab.render_movie" title="qlab.render_movie"><code class="xref py py-func docutils literal notranslate"><span class="pre">render_movie()</span></code></a> function. Arrange the AtomEye window so that the
crack is on the left hand side of the window at the beginning of the simulation
and near the right hand side at the end, then run the command:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">render_movie</span><span class="p">(</span><span class="s1">&#39;movie.mp4&#39;</span><span class="p">)</span>
</pre></div>
</div>
<p>This function renders each frame to a <code class="docutils literal notranslate"><span class="pre">.jpg</span></code> file, before combining the
snapshots with the <a class="reference external" href="http://www.ffmpeg.org/">ffmpeg</a> tool to make a movie like
this one:</p>
<center>
<video width="720" height="360" controls="controls" poster="http://www.jrkermode.co.uk/_movies/classical-111-poster.jpg">
  <source src="http://www.jrkermode.co.uk/_movies/classical-111.mp4"  type='video/mp4' />
  <source src="http://www.jrkermode.co.uk/_movies/classical-111.ogv"  type='video/ogg; codecs="theora, vorbis"'' />
  <source src="http://www.jrkermode.co.uk/_movies/classical-111.webm" type='video/webm; codecs="vp8.0, vorbis"' />
  <p><b>Your browser does not support HTML5 video.
  <a href="http://www.jrkermode.co.uk/_movies/classical-111.mp4">Download</a> the video instead.
  </b></p>
</video>
</center><p>The example movie above makes the ductile nature of the fracture propagation
much clearer. We see local amorphisation, the formation of
strange <em>sp</em><sup>2</sup> tendrils, and temporary crack arrest. Comparing
again with the <a class="reference internal" href="adaptive-qmmm-theory.html#si-tem-images"><span class="std std-ref">experimental TEM images</span></a> makes it clear
that, as a description of fracture in real silicon, the SW potential falls some
way short.</p>
</div>
<div class="section" id="position-of-the-crack-tip">
<span id="position-crack-tip"></span><h3>Position of the crack tip<a class="headerlink" href="#position-of-the-crack-tip" title="Permalink to this headline">¶</a></h3>
<p>The <a class="reference internal" href="../crack.html#quippy.crack.find_crack_tip_stress_field" title="quippy.crack.find_crack_tip_stress_field"><code class="xref py py-func docutils literal notranslate"><span class="pre">find_crack_tip_stress_field()</span></code></a> function works by
fitting per-atom stresses calculated with the SW potential (the
concept of per-atom stresses will be discussed in more detail below)
in the region near the crack tip to the Irwin solution for a singular
crack tip under Mode I loading, which is of the form</p>
<div class="math notranslate nohighlight">
\[\sigma_{ij}(r, \theta) = \frac{K_I}{2\pi r} f_{ij}(\theta)\]</div>
<p>where <span class="math notranslate nohighlight">\(K_I\)</span> is the Mode I stress intensity factor, and the
angular dependence is given by the set of universal functions
<span class="math notranslate nohighlight">\(f_{ij}(\theta)\)</span>.</p>
<p>You can verify this by comparing the position detected by
<a class="reference internal" href="../crack.html#quippy.crack.find_crack_tip_stress_field" title="quippy.crack.find_crack_tip_stress_field"><code class="xref py py-func docutils literal notranslate"><span class="pre">find_crack_tip_stress_field()</span></code></a>,  stored in the
<cite>crack_pos</cite> attribute, with the positions of atoms that visually look
to be near the tip — <cite>right click</cite> on atoms in the AtomEye
viewer window to print information about them, including their
positions.</p>
<p>Compare the automatically detected crack position (printed as the
<cite>crack_pos_x</cite> parameter when you change frames in the AtomEye viewer,
or available via <code class="docutils literal notranslate"><span class="pre">gcat().info['crack_pos_x']</span></code>) with what a visual
inspection of the crack system would tell you. Do you think it’s
accurate enough to use as the basis for selecting a region around the
crack tip to be treated at the QM level?</p>
</div>
<div class="section" id="evolution-of-energy-release-rate-and-crack-position">
<span id="plot-g-and-crack-pos-x"></span><h3>Evolution of energy release rate and crack position<a class="headerlink" href="#evolution-of-energy-release-rate-and-crack-position" title="Permalink to this headline">¶</a></h3>
<p>For <a class="reference internal" href="../io.html#netcdf"><span class="std std-ref">NetCDF</span></a> trajectories,
the <code class="xref py py-attr docutils literal notranslate"><span class="pre">AtomsReaderViewer.reader.netcdf_file</span></code> attribute of the current
viewer object <a class="reference internal" href="../qlab.html#qlab.gcv" title="qlab.gcv"><code class="xref py py-func docutils literal notranslate"><span class="pre">gcv()</span></code></a> provides direct access to the underlying NetCDF
file using the Python <a class="reference external" href="http://code.google.com/p/netcdf4-python/">netCDF4 module</a>:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">traj</span> <span class="o">=</span> <span class="n">gcv</span><span class="p">()</span>
<span class="n">dataset</span> <span class="o">=</span> <span class="n">traj</span><span class="o">.</span><span class="n">reader</span><span class="o">.</span><span class="n">netcdf_file</span>
</pre></div>
</div>
<p>You can list the variables stored in <cite>dataset</cite> with:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="nb">print</span> <span class="n">dataset</span><span class="o">.</span><span class="n">variables</span><span class="o">.</span><span class="n">keys</span><span class="p">()</span>
</pre></div>
</div>
<p>To plot the energy release rate <cite>G</cite> as a function of simulation time,
you could do:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">plot</span><span class="p">(</span><span class="n">dataset</span><span class="o">.</span><span class="n">variables</span><span class="p">[</span><span class="s1">&#39;time&#39;</span><span class="p">],</span> <span class="n">dataset</span><span class="o">.</span><span class="n">variables</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">])</span>
</pre></div>
</div>
<p>You should see that the energy release rate increases at a roughly
constant rate before stopping at constant value when the crack starts
to move (the increase is not linear since is is actually the <cite>strain</cite>
that we increment at a constant rate).</p>
<p>The following plot shows the evolution of <cite>G</cite> (blue) and of the
position of the crack (red; stored as <cite>crack_pos_x</cite>). Note that a
second vertical axis can be produced with the
<a class="reference external" href="https://matplotlib.org/api/_as_gen/matplotlib.pyplot.twinx.html#matplotlib.pyplot.twinx" title="(in Matplotlib v3.0.3)"><code class="xref py py-func docutils literal notranslate"><span class="pre">twinx()</span></code></a> function.</p>
<a class="reference internal image-reference" href="../_images/energy-release-rate-crack-position.png"><img alt="../_images/energy-release-rate-crack-position.png" class="align-center" src="../_images/energy-release-rate-crack-position.png" style="width: 600px;" /></a>
<p>In this case the crack actually arrests for a while at around <span class="math notranslate nohighlight">\(t
= 6\)</span> ps. This is another characteristic feature of non-brittle
fracture, indicating that our simulation is failing to match well
with experiment. According to Griffith’s criterion, fracture should
initiate at <span class="math notranslate nohighlight">\(2\gamma \sim 2.7\)</span> J/m<sup>2</sup>, whereas we
don’t see any motion of the crack tip until <span class="math notranslate nohighlight">\(G \sim 11\)</span> J/m<sup>2</sup>. How much of this difference do you think is due to
the high strain rate and small system used here, and how much to the
choice of interatomic potential? How would you check this?</p>
</div>
<div class="section" id="temperature-and-velocity-analysis">
<span id="plot-temperature"></span><h3>Temperature and velocity analysis<a class="headerlink" href="#temperature-and-velocity-analysis" title="Permalink to this headline">¶</a></h3>
<p>Using the method above, plot the evolution of the temperature during
your simulation. Here is another example plot, with the temperature
shown in blue and the crack position in red.</p>
<a class="reference internal image-reference" href="../_images/temperature-crack-position.png"><img alt="../_images/temperature-crack-position.png" class="align-center" src="../_images/temperature-crack-position.png" style="width: 600px;" /></a>
<p>You will see that lots of heat is produced once the crack starts to
move, indicating that the system is far from equilibrium. This is
another sign that our system is rather small and our strain rate is
rather high. How could this be addressed? Do you think an NVT
simulation would be more realistic? What problems could adding a
thermostat introduce?</p>
<p>If you have time, you could compare how well the atomic velocities
match the expected Maxwell-Boltzmann distribution of atomic
velocities, given by</p>
<div class="math notranslate nohighlight">
\[f(v)\,\mathrm{d}v = 4 \pi \left( \frac{m}{2 \pi k_B T} \right)^{3/2} v^2 \exp \left[ -\frac{mv^2}{2 k_B T} \right] \mathrm{d}v\]</div>
<p>Here’s a Python function which implements this formula:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="k">def</span> <span class="nf">max_bolt</span><span class="p">(</span><span class="n">m</span><span class="p">,</span><span class="n">T</span><span class="p">,</span><span class="n">v</span><span class="p">):</span>
   <span class="s2">&quot;Maxwell-Boltmann distribution of speeds at temperature T for particles of mass m&quot;</span>
   <span class="k">return</span> <span class="mi">4</span><span class="o">*</span><span class="n">pi</span><span class="o">*</span><span class="p">(</span><span class="n">m</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">kB</span><span class="o">*</span><span class="n">T</span><span class="p">))</span><span class="o">**</span><span class="p">(</span><span class="mf">3.0</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">v</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">m</span><span class="o">*</span><span class="n">v</span><span class="o">**</span><span class="mi">2</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">units</span><span class="o">.</span><span class="n">kB</span><span class="o">*</span><span class="n">T</span><span class="p">))</span>
</pre></div>
</div>
<p>We can average the atomic speeds in the last 50 frames in our
trajectory and use the speeds data to produce a histogram:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">m</span> <span class="o">=</span> <span class="n">traj</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">get_masses</span><span class="p">()[</span><span class="mi">0</span><span class="p">]</span>      <span class="c1"># Mass of a Si atom</span>
<span class="n">T</span> <span class="o">=</span> <span class="n">traj</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;temperature&#39;</span><span class="p">]</span>  <span class="c1"># Temperature at end of simulation</span>
<span class="n">v</span> <span class="o">=</span> <span class="n">traj</span><span class="o">.</span><span class="n">reader</span><span class="o">.</span><span class="n">netcdf_file</span><span class="o">.</span><span class="n">variables</span><span class="p">[</span><span class="s1">&#39;momenta&#39;</span><span class="p">][</span><span class="o">-</span><span class="mi">50</span><span class="p">:,:,:]</span><span class="o">/</span><span class="n">m</span> <span class="c1"># Get velocities</span>
<span class="n">s</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">((</span><span class="n">v</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">2</span><span class="p">))</span>      <span class="c1"># Speeds are magnitude of velocities</span>

<span class="n">hist</span><span class="p">(</span><span class="n">s</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">),</span> <span class="mi">20</span><span class="p">,</span> <span class="n">normed</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">alpha</span><span class="o">=</span><span class="mf">0.5</span><span class="p">)</span>  <span class="c1"># Draw a histogram</span>

<span class="n">ss</span> <span class="o">=</span> <span class="n">linspace</span><span class="p">(</span><span class="mf">0.</span><span class="p">,</span> <span class="n">s</span><span class="o">.</span><span class="n">max</span><span class="p">(),</span> <span class="mi">100</span><span class="p">)</span>  <span class="c1"># Compare with Maxwell-Boltzmann distrib</span>
<span class="n">plot</span><span class="p">(</span><span class="n">ss</span><span class="p">,</span> <span class="n">max_bolt</span><span class="p">(</span><span class="n">m</span><span class="p">,</span><span class="n">T</span><span class="p">,</span><span class="n">ss</span><span class="p">),</span> <span class="n">lw</span><span class="o">=</span><span class="mi">2</span><span class="p">)</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="../_images/crack-max-bolt-distrib.png"><img alt="../_images/crack-max-bolt-distrib.png" class="align-center" src="../_images/crack-max-bolt-distrib.png" style="width: 600px;" /></a>
</div>
<div class="section" id="atom-resolved-strain-tensor">
<span id="arsf"></span><h3>Atom-resolved strain tensor<a class="headerlink" href="#atom-resolved-strain-tensor" title="Permalink to this headline">¶</a></h3>
<p>The virial stress expression above is only valid when averaged over
time and space, so this method of calculating per-atom stresses can
lead to unphysical oscillations <a class="reference internal" href="adaptive-qmmm-references.html#zimmerman2004" id="id2">[Zimmerman2004]</a>. One alternative is the
atom-resolved strain tensor, which allows the strain, and hence stress,
fields to be evaluated at the atomistic scale facilitating direct
comparisons with elasticity theory results <a class="reference internal" href="adaptive-qmmm-references.html#moras2010" id="id3">[Moras2010]</a>.</p>
<p>A definition of the atom-resolved strain tensor can be obtained for
all the four-fold coordinated atoms in the tetrahedral structure (all
other atoms are assigned zero strain) by comparing the atomic
positions with the unstrained crystal. The neighbours of each atom are
used to define a local set of cubic axes, and the deformations along
each of these axes are combined into a matrix <span class="math notranslate nohighlight">\(E\)</span> describing the
local deformation:</p>
<div class="math notranslate nohighlight">
\[\begin{split}E = \left(\begin{array}{ccc}
| &amp; | &amp; | \\
\mathbf{e}_{1} &amp; \mathbf{e}_{2} &amp; \mathbf{e}_{3} \\
| &amp; | &amp; |
\end{array}\right)\end{split}\]</div>
<p>where, for example <span class="math notranslate nohighlight">\(\mathbf{e}_{1}\)</span> is the relative deformation
along the first cubic axis.  To compute the local strain of the atom,
we need to separate this deformation into a contribution due to
rotation and one due to strain.  This can be done by finding the polar
decomposition of <span class="math notranslate nohighlight">\(E\)</span>, by writing <span class="math notranslate nohighlight">\(E\)</span> in the form <span class="math notranslate nohighlight">\(E
= SR\)</span> with <span class="math notranslate nohighlight">\(R\)</span> a pure rotation and <span class="math notranslate nohighlight">\(S\)</span> a symmetric matrix.</p>
<p>Diagonalising the product <span class="math notranslate nohighlight">\(EE^T\)</span> allows <span class="math notranslate nohighlight">\(R\)</span> and <span class="math notranslate nohighlight">\(S\)</span>
to be calculated. The strain components <span class="math notranslate nohighlight">\(\epsilon_{xx}\)</span>,
<span class="math notranslate nohighlight">\(\epsilon_{yy}\)</span>, <span class="math notranslate nohighlight">\(\epsilon_{zz}\)</span>, <span class="math notranslate nohighlight">\(\epsilon_{xy}\)</span>,
<span class="math notranslate nohighlight">\(\epsilon_{xz}\)</span> and <span class="math notranslate nohighlight">\(\epsilon_{yz}\)</span> can then be calculated
by rotating <span class="math notranslate nohighlight">\(S\)</span> to align the local cubic axes with the Cartesian
axes:</p>
<div class="math notranslate nohighlight">
\[\begin{split}  R^T S R = I + \epsilon = \left(\begin{array}{ccc}
1 + \epsilon_{xx} &amp; \frac{1}{2}\epsilon_{xy} &amp; \frac{1}{2}\epsilon_{xz} \\
\frac{1}{2}\epsilon_{xy} &amp; 1 + \epsilon_{yy} &amp; \frac{1}{2}\epsilon_{yz} \\
\frac{1}{2}\epsilon_{xz} &amp; \frac{1}{2}\epsilon_{yz} &amp; 1 + \epsilon_{zz}
\end{array}\right).\end{split}\]</div>
<p>Finally if we assume linear elasticity applies, the atomistic stress
can be computed simply as <span class="math notranslate nohighlight">\(\bm\sigma = C \bm\epsilon\)</span> where
<span class="math notranslate nohighlight">\(C\)</span> is the <span class="math notranslate nohighlight">\(6\times6\)</span> matrix of elastic constants.</p>
<p>The <a class="reference internal" href="../elasticity.html#quippy.elasticity.AtomResolvedStressField" title="quippy.elasticity.AtomResolvedStressField"><code class="xref py py-class docutils literal notranslate"><span class="pre">AtomResolvedStressField</span></code></a> class
implements this approach. To use it to calculate the stress in your
<cite>crack_slab</cite> Atoms object, you can use the following code:</p>
<div class="highlight-default notranslate"><div class="highlight"><pre><span></span><span class="n">arsf</span> <span class="o">=</span> <span class="n">AtomResolvedStressField</span><span class="p">(</span><span class="n">bulk</span><span class="o">=</span><span class="n">si_bulk</span><span class="p">)</span>
<span class="n">crack_slab</span><span class="o">.</span><span class="n">set_calculator</span><span class="p">(</span><span class="n">arsf</span><span class="p">)</span>
<span class="n">ar_stress</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">get_stresses</span><span class="p">()</span>
</pre></div>
</div>
<p>Colour your atoms by the <span class="math notranslate nohighlight">\(\sigma_{yy}\)</span> component of the
atom-resolved stress field, and compare with the local virial stress
results. Add the atom resolved <span class="math notranslate nohighlight">\(\sigma_{yy}\)</span> values along
<span class="math notranslate nohighlight">\(y = 0\)</span> to your plot. Do you notice any significant differences?
Repeat the minimisation of the crack slab with a lower value of
<cite>relax_fmax</cite> (e.g. <span class="math notranslate nohighlight">\(1 \times 10^{-3}\)</span> eV/A). Do the stress
components computed using the two methods change much?</p>
<p>When you are ready, proceed to <a class="reference internal" href="adaptive-qmmm-step3.html#step3"><span class="std std-ref">Step 3: LOTF hybrid MD simulation of fracture in Si</span></a>.</p>
</div>
</div>
</div>


           </div>
           
          </div>
          <footer>
  
    <div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
      
        <a href="adaptive-qmmm-step3.html" class="btn btn-neutral float-right" title="Step 3: LOTF hybrid MD simulation of fracture in Si" accesskey="n" rel="next">Next <span class="fa fa-arrow-circle-right"></span></a>
      
      
        <a href="adaptive-qmmm-step1.html" class="btn btn-neutral float-left" title="Step 1: Setup of the Silicon model system" accesskey="p" rel="prev"><span class="fa fa-arrow-circle-left"></span> Previous</a>
      
    </div>
  

  <hr/>

  <div role="contentinfo">
    <p>
        &copy; Copyright 2008-2016, James Kermode

    </p>
  </div>
  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>. 

</footer>

        </div>
      </div>

    </section>

  </div>
  


  <script type="text/javascript">
      jQuery(function () {
          SphinxRtdTheme.Navigation.enable(true);
      });
  </script>

  
  
    
   

</body>
</html>